<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>

<title>Example of a prediction model for significant wave height (SWH) based on sea level pressure (SLP)</title>

<script type="text/javascript">
window.onload = function() {
  var imgs = document.getElementsByTagName('img'), i, img;
  for (i = 0; i < imgs.length; i++) {
    img = imgs[i];
    // center an image if it is the only element of its parent
    if (img.parentElement.childElementCount === 1)
      img.parentElement.style.textAlign = 'center';
  }
};
</script>

<!-- Styles for R syntax highlighter -->
<style type="text/css">
   pre .operator,
   pre .paren {
     color: rgb(104, 118, 135)
   }

   pre .literal {
     color: #990073
   }

   pre .number {
     color: #099;
   }

   pre .comment {
     color: #998;
     font-style: italic
   }

   pre .keyword {
     color: #900;
     font-weight: bold
   }

   pre .identifier {
     color: rgb(0, 0, 0);
   }

   pre .string {
     color: #d14;
   }
</style>

<!-- R syntax highlighter -->
<script type="text/javascript">
var hljs=new function(){function m(p){return p.replace(/&/gm,"&amp;").replace(/</gm,"&lt;")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
hljs.initHighlightingOnLoad();
</script>



<style type="text/css">
body, td {
   font-family: sans-serif;
   background-color: white;
   font-size: 13px;
}

body {
  max-width: 800px;
  margin: auto;
  padding: 1em;
  line-height: 20px;
}

tt, code, pre {
   font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
}

h1 {
   font-size:2.2em;
}

h2 {
   font-size:1.8em;
}

h3 {
   font-size:1.4em;
}

h4 {
   font-size:1.0em;
}

h5 {
   font-size:0.9em;
}

h6 {
   font-size:0.8em;
}

a:visited {
   color: rgb(50%, 0%, 50%);
}

pre, img {
  max-width: 100%;
}
pre {
  overflow-x: auto;
}
pre code {
   display: block; padding: 0.5em;
}

code {
  font-size: 92%;
  border: 1px solid #ccc;
}

code[class] {
  background-color: #F8F8F8;
}

table, td, th {
  border: none;
}

blockquote {
   color:#666666;
   margin:0;
   padding-left: 1em;
   border-left: 0.5em #EEE solid;
}

hr {
   height: 0px;
   border-bottom: none;
   border-top-width: thin;
   border-top-style: dotted;
   border-top-color: #999999;
}

@media print {
   * {
      background: transparent !important;
      color: black !important;
      filter:none !important;
      -ms-filter: none !important;
   }

   body {
      font-size:12pt;
      max-width:100%;
   }

   a, a:visited {
      text-decoration: underline;
   }

   hr {
      visibility: hidden;
      page-break-before: always;
   }

   pre, blockquote {
      padding-right: 1em;
      page-break-inside: avoid;
   }

   tr, img {
      page-break-inside: avoid;
   }

   img {
      max-width: 100% !important;
   }

   @page :left {
      margin: 15mm 20mm 15mm 10mm;
   }

   @page :right {
      margin: 15mm 10mm 15mm 20mm;
   }

   p, h2, h3 {
      orphans: 3; widows: 3;
   }

   h2, h3 {
      page-break-after: avoid;
   }
}
</style>



</head>

<body>
<h2>Example of a prediction model for significant wave height (SWH) based on sea level pressure (SLP)</h2>

<h3>Authors: Hugo L. Hammer, Tor Arne Oigaard, Thordis Thorarinsdottir and Hanne Rognebakke</h3>

<p>Scripts are located in the folder script.
Data example is located in the folder data. </p>

<p>R working directory should be set to the location of the script folder.</p>

<p>To run the analysis, run the script run_models.r as shown below. </p>

<p>Load necessary packages and source functions.</p>

<pre><code class="r">require(glmnet)
require(forecast)
require(moments)
source(&quot;getPreddistr.r&quot;)
source(&quot;BoxCoxLambda.r&quot;)
source(&quot;BoxCoxLambdaKnown.r&quot;)
source(&quot;fourier.r&quot;)
source(&quot;qBoxCox.r&quot;)
source(&quot;pBoxCox.r&quot;)
</code></pre>

<p>Load ERA-Interim data </p>

<pre><code class="r">load(&quot;../data/SWH_SLP_SLP.grad.subset.7.Rdata&quot;)
longitudeSWH
</code></pre>

<pre><code>## [1] -58.50 -57.75 -57.00 -56.25 -55.50 -54.75 -54.00
</code></pre>

<pre><code class="r">latitudeSWH
</code></pre>

<pre><code>## [1] 57.75 57.00 56.25 55.50 54.75 54.00 53.25
</code></pre>

<pre><code class="r">range(years.all)
</code></pre>

<pre><code>## [1] 2006 2015
</code></pre>

<pre><code class="r">summary(SWH)
</code></pre>

<pre><code>##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4726  2.1296  3.0538  3.4142  4.3071 16.5310
</code></pre>

<pre><code class="r">summary(SLP)
</code></pre>

<pre><code>##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   93301  100127  101100  100975  101975  104539
</code></pre>

<pre><code class="r">summary(SLP.grad)
</code></pre>

<pre><code>##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.004531 0.011082 0.017110 0.022871 0.538579
</code></pre>

<p>Set which grid cell to model (out of a total of 7x7 cells).</p>

<pre><code class="r">lonSWH = 4
latSWH = 4
</code></pre>

<p>Define how many missing values is &ldquo;good enough&rdquo;, and define the model to use.</p>

<pre><code class="r">na.thresh = 500
spatial.neighborhoods = 2 
models.fourier = c(&quot;Proposed model, AR order = 5, spatial neigb. = &quot;, spatial.neighborhoods, &quot;\n&quot;)
n.models.fourier = length(models.fourier)
</code></pre>

<p>Create Fourier terms to use in the modeling</p>

<pre><code class="r">m.fourier = 365.25*4 # Number of observations per year. Four observations per day.
K.fourier = 2 # Number of Fourier terms
n.time = length(time.all)
intercept.fourier = fourier(rep(1,n.time), K=K.fourier, m=m.fourier)
colnames(intercept.fourier) = paste(&quot;intercept&quot;, colnames(intercept.fourier), sep = &quot;_&quot;)
</code></pre>

<p>Split that data in a training set and a test set </p>

<pre><code class="r">training.test = list()
training.test[[1]] = which(years.all %in% 2006:2014)
training.test[[2]] = which(years.all %in% 2015)
</code></pre>

<p>Estimate model parameters and obtain predictive distributions</p>

<pre><code class="r">pred.dist  &lt;- getPreddistr(j=lonSWH, k=latSWH, neig=spatial.neighborhoods)
</code></pre>

<pre><code>## Fourier: j = 4 , k = 4 
## Total number of potential predictors: 69 
## Number of predictors selected by LASSO: 52
</code></pre>

<pre><code class="r">pred.mean = pred.dist$pred.mean
pred.sd = pred.dist$pred.sd
pred.lambda = pred.dist$pred.lambda
print(pred.dist$fits)
</code></pre>

<pre><code>## 
## Call:
## lm(formula = SWH.bc.standard.training ~ ., data = predictors.training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76520 -0.10036 -0.00888  0.09102  1.03424 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(&gt;|t|)    
## (Intercept)  8.826e-03  1.912e-03   4.615 3.97e-06 ***
## V1           3.879e-01  2.789e-02  13.908  &lt; 2e-16 ***
## V2          -3.093e-01  2.319e-02 -13.340  &lt; 2e-16 ***
## V4          -2.570e-02  8.679e-03  -2.961 0.003073 ** 
## V5          -1.649e-02  4.442e-03  -3.712 0.000206 ***
## V6           1.227e-01  1.683e-02   7.291 3.26e-13 ***
## V7          -7.061e-06  2.143e-02   0.000 0.999737    
## V8          -1.418e-02  3.901e-02  -0.364 0.716153    
## V9          -6.143e-04  2.416e-03  -0.254 0.799308    
## V10         -5.574e-03  2.476e-03  -2.251 0.024423 *  
## V11         -8.604e-03  2.406e-02  -0.358 0.720678    
## V12          8.637e-02  3.192e-02   2.706 0.006828 ** 
## V20         -3.691e-02  1.218e-02  -3.031 0.002443 ** 
## V22          1.210e-02  1.165e-02   1.039 0.298862    
## V23         -3.883e-03  4.536e-03  -0.856 0.392043    
## V24          3.937e-04  6.237e-03   0.063 0.949668    
## V25          7.161e-03  4.507e-03   1.589 0.112058    
## V26          8.851e-03  6.078e-03   1.456 0.145343    
## V27         -4.226e-03  6.881e-03  -0.614 0.539134    
## V28          5.444e-02  2.339e-02   2.327 0.019956 *  
## V29         -1.269e-02  9.536e-03  -1.331 0.183214    
## V30         -3.398e-02  2.142e-02  -1.587 0.112650    
## V31          1.456e-02  2.289e-03   6.362 2.06e-10 ***
## V32          7.912e-02  2.908e-03  27.213  &lt; 2e-16 ***
## V33         -1.493e-03  2.373e-03  -0.629 0.529159    
## V34         -7.631e-03  2.378e-03  -3.209 0.001336 ** 
## V35         -1.345e-01  3.831e-02  -3.511 0.000447 ***
## V36         -1.921e-02  7.439e-03  -2.582 0.009820 ** 
## V38         -7.006e-02  5.345e-02  -1.311 0.189992    
## V42         -2.294e-02  1.023e-02  -2.242 0.024995 *  
## V43         -1.362e-03  9.077e-03  -0.150 0.880712    
## V44          1.945e-02  1.032e-02   1.884 0.059577 .  
## V45          1.222e+00  8.132e-03 150.238  &lt; 2e-16 ***
## V46          7.960e-03  6.465e-03   1.231 0.218264    
## V47         -1.204e-03  1.149e-02  -0.105 0.916519    
## V48          2.209e-02  9.986e-03   2.212 0.026970 *  
## V49         -1.522e-02  1.137e-02  -1.339 0.180465    
## V50         -4.746e-01  1.339e-02 -35.453  &lt; 2e-16 ***
## V52         -2.856e-02  1.887e-02  -1.514 0.130104    
## V53         -2.621e-02  1.061e-02  -2.470 0.013525 *  
## V54          1.403e-02  1.843e-02   0.761 0.446469    
## V55          1.016e-01  1.397e-02   7.274 3.69e-13 ***
## V56         -2.348e-02  1.371e-02  -1.712 0.086869 .  
## V57          9.316e-03  1.968e-02   0.473 0.635938    
## V59         -1.008e-03  1.365e-02  -0.074 0.941149    
## V60         -4.699e-02  1.303e-02  -3.606 0.000312 ***
## V61          2.440e-02  1.783e-02   1.368 0.171242    
## V62          1.710e-02  1.835e-02   0.932 0.351247    
## V65          2.277e-02  7.117e-03   3.199 0.001382 ** 
## V66         -4.193e-03  1.007e-02  -0.416 0.677244    
## V67         -9.406e-03  1.002e-02  -0.939 0.347790    
## V68          8.730e-03  3.923e-03   2.225 0.026091 *  
## V69         -2.900e-03  5.471e-03  -0.530 0.596090    
## ---
## Signif. codes:  0 &#39;***&#39; 0.001 &#39;**&#39; 0.01 &#39;*&#39; 0.05 &#39;.&#39; 0.1 &#39; &#39; 1
## 
## Residual standard error: 0.1601 on 13095 degrees of freedom
## Multiple R-squared:  0.9745, Adjusted R-squared:  0.9744 
## F-statistic:  9615 on 52 and 13095 DF,  p-value: &lt; 2.2e-16
</code></pre>

<p>Create the vector of observations in the test period</p>

<pre><code class="r">obs  &lt;- SWH[lonSWH, latSWH, training.test[[2]]]
</code></pre>

<p>Calculate PIT values and plot a PIT histogram for the test set </p>

<pre><code class="r">pit  &lt;- pBoxCox(obs, pred.mean, pred.sd, pred.lambda)
hist(pit, freq=FALSE, nclass=10, col=&quot;gray&quot;, xlab=&quot;PIT value&quot;, 
     main=paste(&quot;Lon = &quot;, longitudeSWH[lonSWH], &quot;, Lat = &quot;, latitudeSWH[latSWH]))
abline(a=1, b=0, lty=2)
</code></pre>

<p><img src="" alt="plot of chunk unnamed-chunk-9"/></p>

<p>Calculate the mean absolute error over the test period</p>

<pre><code class="r">pred.median  &lt;- qBoxCox(0.5, pred.mean, pred.sd, pred.lambda)
mae  &lt;- mean(abs(obs - pred.median))
mae
</code></pre>

<pre><code>## [1] 0.2328819
</code></pre>

<p>Plot the prediction and the observation in the first and last 100 time points of the test period</p>

<pre><code class="r">upper  &lt;- qBoxCox(0.95, pred.mean, pred.sd, pred.lambda)
lower  &lt;- qBoxCox(0.05, pred.mean, pred.sd, pred.lambda)

t.period  &lt;- c(1:100)
t.upper  &lt;- upper[t.period]
t.lower  &lt;- lower[t.period]

plot(t.period, obs[t.period], type=&quot;l&quot;,
     xlab=&quot;Time point in test period&quot;, ylab=&quot;SWH&quot;, ylim=c(0,15),
     main=paste(&quot;Lon = &quot;, longitudeSWH[lonSWH], &quot;, Lat = &quot;, latitudeSWH[latSWH]))
polygon(c(t.period, rev(t.period), t.period[1]),
        c(t.lower, rev(t.upper), t.lower[1]), col=&quot;gray90&quot;, border=&quot;NA&quot;)
lines(t.period, obs[t.period])
lines(t.period, pred.median[t.period], col=&quot;gray50&quot;, lty=2)
legend(&quot;topright&quot;, lty=c(1,2,1), lwd=c(1,1,4), col=c(&quot;black&quot;, &quot;gray50&quot;, &quot;gray90&quot;),
       legend=c(&quot;Observation&quot;, &quot;Predictive median&quot;, &quot;90% prediction interval&quot;))
</code></pre>

<p><img src="" alt="plot of chunk unnamed-chunk-11"/></p>

<pre><code class="r">t.period  &lt;- c(1361:1460)
t.upper  &lt;- upper[t.period]
t.lower  &lt;- lower[t.period]

plot(t.period, obs[t.period], type=&quot;l&quot;,
     xlab=&quot;Time point in test period&quot;, ylab=&quot;SWH&quot;, ylim=c(0,15),
     main=paste(&quot;Lon = &quot;, longitudeSWH[lonSWH], &quot;, Lat = &quot;, latitudeSWH[latSWH]))
polygon(c(t.period, rev(t.period), t.period[1]),
        c(t.lower, rev(t.upper), t.lower[1]), col=&quot;gray90&quot;, border=&quot;NA&quot;)
lines(t.period, obs[t.period])
lines(t.period, pred.median[t.period], col=&quot;gray50&quot;, lty=2)
legend(&quot;topright&quot;, lty=c(1,2,1), lwd=c(1,1,4), col=c(&quot;black&quot;, &quot;gray50&quot;, &quot;gray90&quot;),
       legend=c(&quot;Observation&quot;, &quot;Predictive median&quot;, &quot;90% prediction interval&quot;))
</code></pre>

<p><img src="" alt="plot of chunk unnamed-chunk-11"/></p>

<p>Random predictive trajectories for 10 time points</p>

<pre><code class="r">t.ind  &lt;- c(40:49)
random.q  &lt;- array(NA, dim=c(10,10))
for(i in 1:10) random.q[,i]  &lt;- qBoxCox(runif(10), pred.mean[t.ind[i]], pred.sd, pred.lambda)
random.q
</code></pre>

<pre><code>##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
##  [1,] 9.429771  9.610763  8.525877 10.158165 10.332912 10.492195 11.182405
##  [2,] 8.653724  9.316949  9.590308  8.772157  9.875716  9.957759  8.726901
##  [3,] 9.333055  9.577181 10.031716  9.825269  9.581286  9.467095 10.182503
##  [4,] 9.328564 11.045167  9.441738  8.446240  9.158470  7.760592 10.253491
##  [5,] 9.168353  9.088379  9.428345  9.928440  8.635019 10.421799  8.849250
##  [6,] 8.466053 10.536718  9.782850  8.466323  9.471122  9.146203  9.705457
##  [7,] 8.558771  9.293452  8.226525  8.478779  9.976094 10.945288 10.270651
##  [8,] 8.449548  8.853060  8.294006  9.390250  9.696294  8.579420 10.054352
##  [9,] 9.224228 10.098123  9.865008  9.401214  9.284202  9.074849 10.452984
## [10,] 9.039912  8.940767  9.736257  7.923322  9.181679  9.693540 10.787106
##           [,8]     [,9]    [,10]
##  [1,] 8.931922 8.208266 7.856963
##  [2,] 9.374602 8.957247 8.600213
##  [3,] 9.751878 8.597738 7.824361
##  [4,] 8.772154 7.012031 9.086296
##  [5,] 9.264266 8.349947 9.183543
##  [6,] 9.360806 8.352399 7.549240
##  [7,] 8.925110 8.339122 7.631270
##  [8,] 9.452692 7.923413 8.449031
##  [9,] 8.875158 9.278930 8.140958
## [10,] 9.651270 9.047410 8.379112
</code></pre>

<pre><code class="r">plot(t.ind, random.q[1,], type=&quot;l&quot;, col=&quot;gray50&quot;,
     xlab=&quot;Time point in test period&quot;, ylab=&quot;SWH&quot;, ylim=c(5,12),
     main=paste(&quot;Lon = &quot;, longitudeSWH[lonSWH], &quot;, Lat = &quot;, latitudeSWH[latSWH]))
for(i in 2:10) lines(t.ind, random.q[i,], col=&quot;gray50&quot;)
lines(t.ind, obs[t.ind], col=&quot;black&quot;, lwd=2)
</code></pre>

<p><img src="" alt="plot of chunk unnamed-chunk-12"/></p>

<p>Learn correlation from previous timepoints (last 100 time points in test period)</p>

<pre><code class="r">sample.q  &lt;- array(NA, dim=c(10,10))
for(i in 1:10) sample.q[,i]  &lt;- rank(random.q[,i])
sample.q
</code></pre>

<pre><code>##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   10    7    3   10   10    9   10    4    3     4
##  [2,]    4    5    6    5    8    7    1    7    8     8
##  [3,]    9    6   10    8    6    5    5   10    7     3
##  [4,]    8   10    5    2    2    1    6    1    1     9
##  [5,]    6    3    4    9    1    8    2    5    5    10
##  [6,]    2    9    8    3    5    4    3    6    6     1
##  [7,]    3    4    1    4    9   10    7    3    4     2
##  [8,]    1    1    2    6    7    2    4    8    2     7
##  [9,]    7    8    9    7    4    3    8    2   10     5
## [10,]    5    2    7    1    3    6    9    9    9     6
</code></pre>

<pre><code class="r">T  &lt;- length(training.test[[2]])
h.ind  &lt;- training.test[[2]][(T-99):T]
hist.obs  &lt;- t(array(SWH[lonSWH, latSWH, h.ind], dim=c(10,10)))
hist.q  &lt;- array(NA, dim=c(10,10))
for(i in 1:10) hist.q[,i]  &lt;- rank(hist.obs[,i])
hist.q
</code></pre>

<pre><code>##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    5    6    8    8    5    5    8    8    6     8
##  [2,]    8    7    7    7    7    7   10    9    7     4
##  [3,]    2    1    1    1    1    1    1    1    1     2
##  [4,]    3    2    2    2    2    2    2    2    2     3
##  [5,]    4    4    3    3    3    3    4    6    5     9
##  [6,]   10    8    6    6    6    6    5    4    8     6
##  [7,]    7    9   10   10   10   10    9   10    9     7
##  [8,]    6    5    4    4    4    4    3    3    3     1
##  [9,]    1    3    5    5    9    9    6    5    4     5
## [10,]    9   10    9    9    8    8    7    7   10    10
</code></pre>

<pre><code class="r">sort.q  &lt;- random.q
for(i in 1:10) sort.q[,i]  &lt;- sort(random.q[,i])[hist.q[,i]]
plot(t.ind, sort.q[1,], type=&quot;l&quot;, col=&quot;gray50&quot;,
     xlab=&quot;Time point in test period&quot;, ylab=&quot;SWH&quot;, ylim=c(5,12),
     main=paste(&quot;Lon = &quot;, longitudeSWH[lonSWH], &quot;, Lat = &quot;, latitudeSWH[latSWH]))
for(i in 2:10) lines(t.ind, sort.q[i,], col=&quot;gray50&quot;)
lines(t.ind, obs[t.ind], col=&quot;black&quot;, lwd=2)
</code></pre>

<p><img src="" alt="plot of chunk unnamed-chunk-13"/></p>

</body>

</html>
